Exploratory analysis

rh

rh

main <- main[!is.na(main$icustay_id), ] # remove 1071 w/o ICU admissions
main <- main[!is.na(main$rh), ] # remove 902 w/o two or more glucose
main$gender <- ifelse(main$gender=="F", 0, 1)

for (i in 1:nrow(main)) {
  if (main$ins[i]==0 & is.na(main$ins_iv[i])) {
    main$ins_iv[i] <- 0
    main$ins_sc[i] <- 0
  } else if (main$ins[i]==1 & is.na(main$ins_iv[i])) {
    main$ins_iv[i] <- 0
    main$ins_sc[i] <- 1 # for those in chartevent but no mention about route -> subcutaneous
  }
}

main2 <- convert.bin.fac(main)

main2$age_cat <- cut2(main2$age, g=4)
main2$sofa_cat <- cut2(main2$sofa_1, g=4)
main2$rh_cat <- cut2(main2$rh_freq, c(0, 1, 3, 5))
main2$rrh_cat <- cut2(main2$rrh_freq, c(0, 1, 3, 5))
#main2$rh_cat <- cut2(main2$rh_freq, g=3)
#main2$rrh_cat <- cut2(main2$rrh_freq, g=3)
table(main2$rh_cat)
## 
##       0 [ 1, 3) [ 3, 5) [ 5,19] 
##   39377    8432     828     186
table(main2$rrh_cat)
## 
##       0 [ 1, 3) [ 3, 5) [ 5,11] 
##   42338    6047     375      63
summary(main2)
##     hadm_id         icustay_id       subject_id         age      
##  Min.   :100001   Min.   :200001   Min.   :    3   Min.   :18.0  
##  1st Qu.:124924   1st Qu.:225040   1st Qu.:13258   1st Qu.:53.0  
##  Median :150126   Median :250114   Median :26672   Median :66.0  
##  Mean   :150015   Mean   :250068   Mean   :36782   Mean   :64.1  
##  3rd Qu.:175086   3rd Qu.:275104   3rd Qu.:60144   3rd Qu.:78.0  
##  Max.   :199999   Max.   :299999   Max.   :99999   Max.   :91.4  
##                                                                  
##  gender    dm           type_dm       a1c_pre_icu        sofa_1      
##  0:21365   0:30755   Min.   :1.00    Min.   : 3.8    Min.   : 0.000  
##  1:27458   1:18068   1st Qu.:2.00    1st Qu.: 5.6    1st Qu.: 0.000  
##                      Median :2.00    Median : 6.0    Median : 2.000  
##                      Mean   :1.88    Mean   : 6.5    Mean   : 2.924  
##                      3rd Qu.:2.00    3rd Qu.: 6.8    3rd Qu.: 5.000  
##                      Max.   :2.00    Max.   :19.0    Max.   :15.000  
##                      NA's   :35123   NA's   :43425                   
##  elixhauser_hospital elixhauser_28day  survival_day       icu_los        
##  Min.   :-25.000     Min.   :-20.00   Min.   :  -3.0   Min.   :  0.0001  
##  1st Qu.: -1.000     1st Qu.:  0.00   1st Qu.:  17.0   1st Qu.:  1.2125  
##  Median :  1.000     Median :  3.00   Median : 145.0   Median :  2.1334  
##  Mean   :  2.486     Mean   :  3.75   Mean   : 512.7   Mean   :  4.0924  
##  3rd Qu.:  6.000     3rd Qu.:  8.00   3rd Qu.: 715.0   3rd Qu.:  4.1369  
##  Max.   : 33.000     Max.   : 38.00   Max.   :4169.0   Max.   :153.9280  
##  NA's   :9           NA's   :9        NA's   :27348    NA's   :3         
##     hosp_los        mort28    mort90    mort365    admit_time       
##  Min.   : -0.8333   0:42244   0:39391   0:35230   Length:48823      
##  1st Qu.:  4.0833   1: 6579   1: 9432   1:13593   Class :character  
##  Median :  6.9583                                 Mode  :character  
##  Mean   : 10.0637                                                   
##  3rd Qu.: 12.0000                                                   
##  Max.   :294.6250                                                   
##                                                                     
##   disch_time         icu_intime        icu_outtime           lac_mean     
##  Length:48823       Length:48823       Length:48823       Min.   : 0.300  
##  Class :character   Class :character   Class :character   1st Qu.: 1.300  
##  Mode  :character   Mode  :character   Mode  :character   Median : 1.800  
##                                                           Mean   : 2.202  
##                                                           3rd Qu.: 2.500  
##                                                           Max.   :24.800  
##                                                           NA's   :15323   
##     lac_max          lac_min           lac_sd       ins       tpn      
##  Min.   : 0.300   Min.   : 0.000   Min.   : 0.000   0:24253   0:46393  
##  1st Qu.: 1.600   1st Qu.: 0.900   1st Qu.: 0.338   1:24570   1: 2430  
##  Median : 2.300   Median : 1.200   Median : 0.636                      
##  Mean   : 3.192   Mean   : 1.518   Mean   : 0.952                      
##  3rd Qu.: 3.600   3rd Qu.: 1.700   3rd Qu.: 1.176                      
##  Max.   :36.000   Max.   :24.100   Max.   :19.021                      
##  NA's   :15323    NA's   :15323    NA's   :24580                       
##  vasop     ventilat     glu_mean        glu_max          glu_min       
##  0:46506   0:24630   Min.   : 29.5   Min.   :  40.0   Min.   :-251.00  
##  1: 2317   1:24193   1st Qu.:109.9   1st Qu.: 143.0   1st Qu.:  75.00  
##                      Median :124.3   Median : 178.0   Median :  87.00  
##                      Mean   :133.4   Mean   : 214.7   Mean   :  88.38  
##                      3rd Qu.:144.8   3rd Qu.: 235.0   3rd Qu.:  99.00  
##                      Max.   :745.1   Max.   :3565.0   Max.   : 621.00  
##                                                                        
##  hypoglycemia hypoglycemia_freq hyperglycemia hyperglycemia_freq rh       
##  0:39359      Min.   : 0.0000   0:25160       Min.   : 0.000     0:39377  
##  1: 9464      1st Qu.: 0.0000   1:23663       1st Qu.: 0.000     1: 9446  
##               Median : 0.0000                 Median : 0.000              
##               Mean   : 0.4173                 Mean   : 2.002              
##               3rd Qu.: 0.0000                 3rd Qu.: 2.000              
##               Max.   :37.0000                 Max.   :97.000              
##                                                                           
##     rh_freq        rrh          rrh_freq       hypo_after_rh
##  Min.   : 0.0000   0:42338   Min.   : 0.0000   0:48823      
##  1st Qu.: 0.0000   1: 6485   1st Qu.: 0.0000                
##  Median : 0.0000             Median : 0.0000                
##  Mean   : 0.2926             Mean   : 0.1791                
##  3rd Qu.: 0.0000             3rd Qu.: 0.0000                
##  Max.   :19.0000             Max.   :11.0000                
##                                                             
##  hypo_after_rrh rh_after_hyper    ins_mean          ins_iv    ins_sc   
##  0:48823        0:48647        Min.   :-10947.000   0:31473   0:33681  
##                 1:  176        1st Qu.:     0.000   1:17350   1:15142  
##                                Median :     2.333                      
##                                Mean   :    -8.787                      
##                                3rd Qu.:     4.691                      
##                                Max.   :   500.000                      
##                                NA's   :24253                           
##       age_cat        sofa_cat         rh_cat         rrh_cat     
##  [18,54.0):12725   0     :14259    0     :39377    0     :42338  
##  [54,67.0):12524   [1, 3):10961   [ 1, 3): 8432   [ 1, 3): 6047  
##  [67,79.0):12086   [3, 6):13884   [ 3, 5):  828   [ 3, 5):  375  
##  [79,91.4]:11488   [6,15]: 9719   [ 5,19]:  186   [ 5,11]:   63  
##                                                                  
##                                                                  
## 
#tapply(main2$dm, main2$mort28, summary)

# CreateTableOne(vars=c("age", "gender", "dm", "type_dm", "sofa_1", "elixhauser_hospital", "ins_mean", "tpn", "vasop", "ventilat", "mort28"),
#                strata="mort28",
#                data=main2, test=TRUE)

Univariate analysis

Diabetes

  1. Demonstrate the baseline covariates and outcomes of our cohort, separated by diabetes (-/+). This might be table 1.
  2. Compute the OR of 28/90/365-day mortality between diabetes(-/+).

    • Conclusion: diabetes group has higher mortality regardless other confounders. This might be figure 2a
    • Diabetes is associated with increasing mortality significantly in 90- and 365-day but not 28-day.
CreateTableOne(vars=c("age", "gender", "a1c_pre_icu", "sofa_1", "elixhauser_hospital", 
                      "hypoglycemia", "hypoglycemia_freq", "hyperglycemia", "hyperglycemia_freq", 
                      "rh", "rh_freq", "rrh", "rrh_freq",
                      "ins", "ins_iv", "ins_sc", "tpn", "vasop", "ventilat", 
                      "mort28", "mort90", "mort365", 
                      "survival_day", "icu_los", "hosp_los"),
               strata=c("dm"),
               data=main2, test=TRUE) %>% print(
  printToggle      = FALSE,
  showAllLevels    = TRUE,
  cramVars         = "kon"
  ) %>% {data.frame(
  variable_name    = gsub(" ", "&nbsp;", rownames(.), fixed = TRUE), ., 
  row.names        = NULL, 
  check.names      = FALSE, 
  stringsAsFactors = FALSE)} %>% knitr::kable()
variable_name level 0 1 p test
n 30755 18068
age (mean (sd)) 62.93 (18.43) 66.09 (14.94) <0.001
gender (%) 0 13408 (43.6) 7957 (44.0) 0.346
1 17347 (56.4) 10111 (56.0)
a1c_pre_icu (mean (sd)) 5.71 (0.41) 7.30 (1.85) <0.001
sofa_1 (mean (sd)) 2.95 (2.83) 2.88 (2.85) 0.014
elixhauser_hospital (mean (sd)) 2.63 (6.27) 2.23 (6.49) <0.001
hypoglycemia (%) 0 26536 (86.3) 12823 (71.0) <0.001
1 4219 (13.7) 5245 (29.0)
hypoglycemia_freq (mean (sd)) 0.28 (1.04) 0.65 (1.57) <0.001
hyperglycemia (%) 0 20125 (65.4) 5035 (27.9) <0.001
1 10630 (34.6) 13033 (72.1)
hyperglycemia_freq (mean (sd)) 0.96 (2.41) 3.77 (5.36) <0.001
rh (%) 0 25926 (84.3) 13451 (74.4) <0.001
1 4829 (15.7) 4617 (25.6)
rh_freq (mean (sd)) 0.22 (0.63) 0.42 (0.94) <0.001
rrh (%) 0 27356 (88.9) 14982 (82.9) <0.001
1 3399 (11.1) 3086 (17.1)
rrh_freq (mean (sd)) 0.14 (0.47) 0.24 (0.63) <0.001
ins (%) 0 18079 (58.8) 6174 (34.2) <0.001
1 12676 (41.2) 11894 (65.8)
ins_iv (%) 0 21822 (71.0) 9651 (53.4) <0.001
1 8933 (29.0) 8417 (46.6)
ins_sc (%) 0 23467 (76.3) 10214 (56.5) <0.001
1 7288 (23.7) 7854 (43.5)
tpn (%) 0 29168 (94.8) 17225 (95.3) 0.016
1 1587 ( 5.2) 843 ( 4.7)
vasop (%) 0 29397 (95.6) 17109 (94.7) <0.001
1 1358 ( 4.4) 959 ( 5.3)
ventilat (%) 0 15453 (50.2) 9177 (50.8) 0.248
1 15302 (49.8) 8891 (49.2)
mort28 (%) 0 26632 (86.6) 15612 (86.4) 0.568
1 4123 (13.4) 2456 (13.6)
mort90 (%) 0 24984 (81.2) 14407 (79.7) <0.001
1 5771 (18.8) 3661 (20.3)
mort365 (%) 0 22625 (73.6) 12605 (69.8) <0.001
1 8130 (26.4) 5463 (30.2)
survival_day (mean (sd)) 501.58 (751.67) 528.60 (743.80) 0.009
icu_los (mean (sd)) 3.95 (5.83) 4.33 (6.26) <0.001
hosp_los (mean (sd)) 9.78 (10.83) 10.54 (10.83) <0.001
# prop.table(table(main2$dm, main2$mort28, dnn=c("DM", "28-day Mortality")), 1)
# prop.table(table(main2$dm, main2$mort90, dnn=c("DM", "90-day Mortality")), 1)
# prop.table(table(main2$dm, main2$mort365, dnn=c("DM", "365-day Mortality")), 1)

CreateTableOne(vars=c("mort28", "mort90", "mort365"),
               strata=c("dm"),
               data=main2, test=TRUE)
##                  Stratified by dm
##                   0             1             p      test
##   n               30755         18068                    
##   mort28 = 1 (%)   4123 (13.4)   2456 (13.6)   0.568     
##   mort90 = 1 (%)   5771 (18.8)   3661 (20.3)  <0.001     
##   mort365 = 1 (%)  8130 (26.4)   5463 (30.2)  <0.001
# par(mfrow=c(1, 3))
# plot(main2$dm, main2$mort28, xlab="DM", ylab="28-day Mortality")
# plot(main2$dm, main2$mort90, xlab="DM", ylab="90-day Mortality")
# plot(main2$dm, main2$mort365, xlab="DM", ylab="365-day Mortality")
# 
# p1 <- plot_prop_by_level(main2, "dm", "mort28")
# p2 <- plot_prop_by_level(main2, "dm", "mort90")
# p3 <- plot_prop_by_level(main2, "dm", "mort365")
# grid.arrange(p1, p2, p3, ncol=3)

o1 <- plot_OR_by_level(main2, "dm", "mort28", include.ref.group.effect=TRUE)
o2 <- plot_OR_by_level(main2, "dm", "mort90", include.ref.group.effect=TRUE)
o3 <- plot_OR_by_level(main2, "dm", "mort365", include.ref.group.effect=TRUE)
grid.arrange(o1, o2, o3, ncol=3)

# t.test(main2[main2$dm==0, ], main2[main2$dm==1, ],
#        alternative="two.sided", var.equal=FALSE, paired=FALSE)
summary(fit.glm <- glm(mort28 ~ dm, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort28 ~ dm, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.5406  -0.5406  -0.5365  -0.5365   2.0047  
## 
## Coefficients:
##             Estimate Std. Error  z value Pr(>|z|)    
## (Intercept) -1.86553    0.01674 -111.469   <2e-16 ***
## dm1          0.01603    0.02741    0.585    0.559    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 38602  on 48822  degrees of freedom
## Residual deviance: 38601  on 48821  degrees of freedom
## AIC: 38605
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                 2.5 %    97.5 %
## (Intercept) 0.1497982 0.1599552
## dm1         0.9629177 1.0721463
sjp.glm(fit.glm)
## Waiting for profiling to be done...

summary(fit.glm <- glm(mort90 ~ dm, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort90 ~ dm, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.6730  -0.6730  -0.6447  -0.6447   1.8293  
## 
## Coefficients:
##             Estimate Std. Error  z value Pr(>|z|)    
## (Intercept) -1.46539    0.01460 -100.335  < 2e-16 ***
## dm1          0.09541    0.02358    4.047 5.19e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 47926  on 48822  degrees of freedom
## Residual deviance: 47910  on 48821  degrees of freedom
## AIC: 47914
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                 2.5 %    97.5 %
## (Intercept) 0.2244502 0.2376753
## dm1         1.0503819 1.1520875
sjp.glm(fit.glm)
## Waiting for profiling to be done...

summary(fit.glm <- glm(mort365 ~ dm, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort365 ~ dm, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.8486  -0.8486  -0.7836   1.5467   1.6313  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.02349    0.01293 -79.153   <2e-16 ***
## dm1          0.18740    0.02073   9.042   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 57753  on 48822  degrees of freedom
## Residual deviance: 57671  on 48821  degrees of freedom
## AIC: 57675
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                 2.5 %    97.5 %
## (Intercept) 0.3503268 0.3685418
## dm1         1.1580697 1.2560875
sjp.glm(fit.glm)
## Waiting for profiling to be done...
## Warning: Removed 1 rows containing missing values (geom_hline).

Relative hypoglycemia

  1. We stratified the admissions into five groups based on the frequency of RH.
  2. Demonstrate the covariates and outcomes of cohort, separated by the frequency of RH (into 4 groups)
  3. Compute the OR of 28/90/365-day mortality between different RH groups.

    • Conclusion: The OR fluctuates. Without RH, the diabetic admissions have higher mortality of 28-, 90- and 365-day. However, it has a trend that the non-diabetic admissions with the event of RH associates with higher mortality (figure 2b). (in discussion) By comarping the odds ratio between the diabetic and non-diabetic group, diabetes becomes a protective factor when the event of RH happens, especially the multiple episodes of RH.
    • RH is associated with increasing mortality in 90-day.
CreateTableOne(vars=c("age", "gender", "a1c_pre_icu", "sofa_1", "elixhauser_hospital", 
                      "hypoglycemia", "hypoglycemia_freq", "hyperglycemia", "hyperglycemia_freq", 
                      "ins", "ins_iv", "ins_sc", "tpn", "vasop", "ventilat", 
                      "mort28", "mort90", "mort365", 
                      "survival_day", "icu_los", "hosp_los"),
               strata=c("rh_cat"),
               data=main2, test=TRUE) %>% print(
  printToggle      = FALSE,
  showAllLevels    = TRUE,
  cramVars         = "kon"
  ) %>% {data.frame(
  variable_name    = gsub(" ", "&nbsp;", rownames(.), fixed = TRUE), ., 
  row.names        = NULL, 
  check.names      = FALSE, 
  stringsAsFactors = FALSE)} %>% knitr::kable()
variable_name level 0 [ 1, 3) [ 3, 5) [ 5,19] p test
n 39377 8432 828 186
age (mean (sd)) 64.28 (17.46) 63.48 (16.54) 62.39 (16.62) 61.89 (15.80) <0.001
gender (%) 0 17259 (43.8) 3610 (42.8) 404 (48.8) 92 ( 49.5) 0.003
1 22118 (56.2) 4822 (57.2) 424 (51.2) 94 ( 50.5)
a1c_pre_icu (mean (sd)) 6.36 (1.34) 6.67 (1.74) 7.40 (2.53) 6.61 (1.19) <0.001
sofa_1 (mean (sd)) 2.57 (2.67) 4.36 (2.97) 4.70 (3.48) 5.10 (3.37) <0.001
elixhauser_hospital (mean (sd)) 2.51 (6.39) 2.33 (6.28) 2.60 (5.84) 4.46 (5.28) <0.001
hypoglycemia (%) 0 33211 (84.3) 5769 (68.4) 358 (43.2) 21 ( 11.3) <0.001
1 6166 (15.7) 2663 (31.6) 470 (56.8) 165 ( 88.7)
hypoglycemia_freq (mean (sd)) 0.30 (0.94) 0.72 (1.69) 1.74 (2.70) 5.89 (5.89) <0.001
hyperglycemia (%) 0 23626 (60.0) 1512 (17.9) 22 ( 2.7) 0 ( 0.0) <0.001
1 15751 (40.0) 6920 (82.1) 806 (97.3) 186 (100.0)
hyperglycemia_freq (mean (sd)) 1.32 (2.79) 4.06 (5.16) 9.25 (8.38) 20.10 (14.85) <0.001
ins (%) 0 21858 (55.5) 2310 (27.4) 84 (10.1) 1 ( 0.5) <0.001
1 17519 (44.5) 6122 (72.6) 744 (89.9) 185 ( 99.5)
ins_iv (%) 0 28217 (71.7) 3143 (37.3) 110 (13.3) 3 ( 1.6) <0.001
1 11160 (28.3) 5289 (62.7) 718 (86.7) 183 ( 98.4)
ins_sc (%) 0 28026 (71.2) 5020 (59.5) 509 (61.5) 126 ( 67.7) <0.001
1 11351 (28.8) 3412 (40.5) 319 (38.5) 60 ( 32.3)
tpn (%) 0 38017 (96.5) 7640 (90.6) 631 (76.2) 105 ( 56.5) <0.001
1 1360 ( 3.5) 792 ( 9.4) 197 (23.8) 81 ( 43.5)
vasop (%) 0 38084 (96.7) 7636 (90.6) 675 (81.5) 111 ( 59.7) <0.001
1 1293 ( 3.3) 796 ( 9.4) 153 (18.5) 75 ( 40.3)
ventilat (%) 0 22458 (57.0) 1986 (23.6) 170 (20.5) 16 ( 8.6) <0.001
1 16919 (43.0) 6446 (76.4) 658 (79.5) 170 ( 91.4)
mort28 (%) 0 34120 (86.6) 7243 (85.9) 728 (87.9) 153 ( 82.3) 0.054
1 5257 (13.4) 1189 (14.1) 100 (12.1) 33 ( 17.7)
mort90 (%) 0 31853 (80.9) 6784 (80.5) 641 (77.4) 113 ( 60.8) <0.001
1 7524 (19.1) 1648 (19.5) 187 (22.6) 73 ( 39.2)
mort365 (%) 0 28354 (72.0) 6220 (73.8) 571 (69.0) 85 ( 45.7) <0.001
1 11023 (28.0) 2212 (26.2) 257 (31.0) 101 ( 54.3)
survival_day (mean (sd)) 497.97 (727.62) 571.63 (820.11) 664.10 (918.84) 368.06 (685.60) <0.001
icu_los (mean (sd)) 3.49 (4.63) 5.77 (7.82) 11.52 (13.88) 22.84 (21.81) <0.001
hosp_los (mean (sd)) 8.96 (9.05) 13.47 (13.48) 21.41 (22.05) 38.62 (30.27) <0.001
# prop.table(table(main2$rh, main2$mort28, dnn=c("RH", "28-day Mortality")), 1)
# prop.table(table(main2$rh, main2$mort90, dnn=c("RH", "90-day Mortality")), 1)
# prop.table(table(main2$rh, main2$mort365, dnn=c("RH", "365-day Mortality")), 1)

CreateTableOne(vars=c("mort28", "mort90", "mort365"),
               strata=c("rh_cat"),
               data=main2, test=TRUE)
##                  Stratified by rh_cat
##                    0            [ 1, 3)      [ 3, 5)     [ 5,19]    
##   n               39377         8432         828         186        
##   mort28 = 1 (%)   5257 (13.4)  1189 (14.1)  100 (12.1)   33 (17.7) 
##   mort90 = 1 (%)   7524 (19.1)  1648 (19.5)  187 (22.6)   73 (39.2) 
##   mort365 = 1 (%) 11023 (28.0)  2212 (26.2)  257 (31.0)  101 (54.3) 
##                  Stratified by rh_cat
##                   p      test
##   n                          
##   mort28 = 1 (%)   0.054     
##   mort90 = 1 (%)  <0.001     
##   mort365 = 1 (%) <0.001
# par(mfrow=c(1, 3))
# plot(main2$rh_cat, main2$mort28, xlab="RH", ylab="28-day Mortality")
# plot(main2$rh_cat, main2$mort90, xlab="RH", ylab="90-day Mortality")
# plot(main2$rh_cat, main2$mort365, xlab="RH", ylab="365-day Mortality")
# 
# p1 <- plot_prop_by_level(main2, "rh_cat", "mort28", factor.var2="dm")
# p2 <- plot_prop_by_level(main2, "rh_cat", "mort90", factor.var2="dm")
# p3 <- plot_prop_by_level(main2, "rh_cat", "mort365", factor.var2="dm")
# grid.arrange(p1, p2, p3, ncol=3)

o1 <- plot_OR_by_level(main2, "rh_cat", "mort28", include.ref.group.effect=TRUE)
o2 <- plot_OR_by_level(main2, "rh_cat", "mort90", include.ref.group.effect=TRUE)
o3 <- plot_OR_by_level(main2, "rh_cat", "mort365", include.ref.group.effect=TRUE)
grid.arrange(o1, o2, o3, ncol=3)

summary(fit.glm <- glm(mort28 ~ rh, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort28 ~ rh, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.5491  -0.5353  -0.5353  -0.5353   2.0068  
## 
## Coefficients:
##             Estimate Std. Error  z value Pr(>|z|)    
## (Intercept) -1.87032    0.01482 -126.232   <2e-16 ***
## rh1          0.05465    0.03315    1.648   0.0993 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 38602  on 48822  degrees of freedom
## Residual deviance: 38599  on 48821  degrees of freedom
## AIC: 38603
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                 2.5 %    97.5 %
## (Intercept) 0.1496484 0.1585975
## rh1         0.9894168 1.1267341
sjp.glm(fit.glm)
## Waiting for profiling to be done...

summary(fit.glm <- glm(mort90 ~ rh, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort90 ~ rh, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.6718  -0.6512  -0.6512  -0.6512   1.8194  
## 
## Coefficients:
##             Estimate Std. Error  z value Pr(>|z|)    
## (Intercept) -1.44303    0.01282 -112.578   <2e-16 ***
## rh1          0.06913    0.02865    2.413   0.0158 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 47926  on 48822  degrees of freedom
## Residual deviance: 47920  on 48821  degrees of freedom
## AIC: 47924
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                 2.5 %    97.5 %
## (Intercept) 0.2303346 0.2422039
## rh1         1.0128623 1.1332717
sjp.glm(fit.glm)
## Waiting for profiling to be done...

summary(fit.glm <- glm(mort365 ~ rh, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort365 ~ rh, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.8105  -0.8105  -0.8105   1.5957   1.6135  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.94478    0.01122 -84.172   <2e-16 ***
## rh1         -0.03935    0.02570  -1.531    0.126    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 57753  on 48822  degrees of freedom
## Residual deviance: 57750  on 48821  degrees of freedom
## AIC: 57754
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                 2.5 %    97.5 %
## (Intercept) 0.3802907 0.3973968
## rh1         0.9140765 1.0109672
sjp.glm(fit.glm)
## Waiting for profiling to be done...

summary(fit.glm <- glm(mort28 ~ rh_cat, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort28 ~ rh_cat, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.6250  -0.5353  -0.5353  -0.5353   2.0561  
## 
## Coefficients:
##               Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)   -1.87032    0.01482 -126.232   <2e-16 ***
## rh_cat[ 1, 3)  0.06340    0.03462    1.831   0.0671 .  
## rh_cat[ 3, 5) -0.11481    0.10767   -1.066   0.2863    
## rh_cat[ 5,19]  0.33639    0.19251    1.747   0.0806 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 38602  on 48822  degrees of freedom
## Residual deviance: 38594  on 48819  degrees of freedom
## AIC: 38602
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                   2.5 %    97.5 %
## (Intercept)   0.1496484 0.1585975
## rh_cat[ 1, 3) 0.9952003 1.1398688
## rh_cat[ 3, 5) 0.7178496 1.0952928
## rh_cat[ 5,19] 0.9442578 2.0132919
sjp.glm(fit.glm)
## Waiting for profiling to be done...

summary(fit.glm <- glm(mort90 ~ rh_cat, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort90 ~ rh_cat, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.9984  -0.6512  -0.6512  -0.6512   1.8194  
## 
## Coefficients:
##               Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)   -1.44303    0.01282 -112.578  < 2e-16 ***
## rh_cat[ 1, 3)  0.02803    0.03031    0.925   0.3550    
## rh_cat[ 3, 5)  0.21111    0.08409    2.510   0.0121 *  
## rh_cat[ 5,19]  1.00611    0.15071    6.676 2.46e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 47926  on 48822  degrees of freedom
## Residual deviance: 47880  on 48819  degrees of freedom
## AIC: 47888
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                   2.5 %    97.5 %
## (Intercept)   0.2303346 0.2422039
## rh_cat[ 1, 3) 0.9688894 1.0911177
## rh_cat[ 3, 5) 1.0447940 1.4530446
## rh_cat[ 5,19] 2.0280928 3.6653657
sjp.glm(fit.glm)
## Waiting for profiling to be done...

summary(fit.glm <- glm(mort365 ~ rh_cat, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort365 ~ rh_cat, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.2515  -0.8105  -0.8105   1.5957   1.6359  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -0.94478    0.01122 -84.172  < 2e-16 ***
## rh_cat[ 1, 3) -0.08909    0.02718  -3.278  0.00105 ** 
## rh_cat[ 3, 5)  0.14647    0.07595   1.929  0.05379 .  
## rh_cat[ 5,19]  1.11725    0.14762   7.568 3.78e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 57753  on 48822  degrees of freedom
## Residual deviance: 57680  on 48819  degrees of freedom
## AIC: 57688
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                   2.5 %    97.5 %
## (Intercept)   0.3802907 0.3973968
## rh_cat[ 1, 3) 0.8671747 0.9646785
## rh_cat[ 3, 5) 0.9962092 1.3418586
## rh_cat[ 5,19] 2.2901436 4.0888922
sjp.glm(fit.glm)
## Waiting for profiling to be done...

Rapid relative hypoglycemia

  1. Demonstrate the covariates and outcomes of cohort, separated by the frequency of RRH (into 4 groups)
  2. Compute the OR of 28/90/365-day mortality between different RRH groups.

    • Conclusion: The OR also fluctuates. The same observation can be found in the episode of RRH (figure 2c).
CreateTableOne(vars=c("age", "gender", "a1c_pre_icu", "sofa_1", "elixhauser_hospital", 
                      "hypoglycemia", "hypoglycemia_freq", "hyperglycemia", "hyperglycemia_freq", 
                      "ins", "ins_iv", "ins_sc", "tpn", "vasop", "ventilat", 
                      "mort28", "mort90", "mort365", 
                      "survival_day", "icu_los", "hosp_los"),
               strata=c("rrh_cat"),
               data=main2, test=TRUE) %>% print(
  printToggle      = FALSE,
  showAllLevels    = TRUE,
  cramVars         = "kon"
  ) %>% {data.frame(
  variable_name    = gsub(" ", "&nbsp;", rownames(.), fixed = TRUE), ., 
  row.names        = NULL, 
  check.names      = FALSE, 
  stringsAsFactors = FALSE)} %>% knitr::kable()
variable_name level 0 [ 1, 3) [ 3, 5) [ 5,11] p test
n 42338 6047 375 63
age (mean (sd)) 64.29 (17.40) 62.90 (16.47) 62.54 (16.26) 62.11 (16.77) <0.001
gender (%) 0 18538 (43.8) 2605 (43.1) 191 (50.9) 31 ( 49.2) 0.021
1 23800 (56.2) 3442 (56.9) 184 (49.1) 32 ( 50.8)
a1c_pre_icu (mean (sd)) 6.42 (1.41) 6.65 (1.75) 7.53 (2.80) 6.72 (1.42) <0.001
sofa_1 (mean (sd)) 2.67 (2.73) 4.55 (2.94) 4.76 (3.40) 5.51 (3.45) <0.001
elixhauser_hospital (mean (sd)) 2.56 (6.40) 1.92 (6.03) 2.47 (5.80) 5.13 (5.19) <0.001
hypoglycemia (%) 0 35198 (83.1) 4019 (66.5) 137 (36.5) 5 ( 7.9) <0.001
1 7140 (16.9) 2028 (33.5) 238 (63.5) 58 ( 92.1)
hypoglycemia_freq (mean (sd)) 0.33 (1.03) 0.82 (1.88) 2.27 (3.40) 7.65 (6.91) <0.001
hyperglycemia (%) 0 24134 (57.0) 1016 (16.8) 10 ( 2.7) 0 ( 0.0) <0.001
1 18204 (43.0) 5031 (83.2) 365 (97.3) 63 (100.0)
hyperglycemia_freq (mean (sd)) 1.56 (3.24) 4.40 (5.70) 10.37 (10.62) 20.11 (15.32) <0.001
ins (%) 0 22767 (53.8) 1451 (24.0) 34 ( 9.1) 1 ( 1.6) <0.001
1 19571 (46.2) 4596 (76.0) 341 (90.9) 62 ( 98.4)
ins_iv (%) 0 29486 (69.6) 1945 (32.2) 40 (10.7) 2 ( 3.2) <0.001
1 12852 (30.4) 4102 (67.8) 335 (89.3) 61 ( 96.8)
ins_sc (%) 0 29794 (70.4) 3598 (59.5) 246 (65.6) 43 ( 68.3) <0.001
1 12544 (29.6) 2449 (40.5) 129 (34.4) 20 ( 31.7)
tpn (%) 0 40677 (96.1) 5401 (89.3) 280 (74.7) 35 ( 55.6) <0.001
1 1661 ( 3.9) 646 (10.7) 95 (25.3) 28 ( 44.4)
vasop (%) 0 40734 (96.2) 5438 (89.9) 303 (80.8) 31 ( 49.2) <0.001
1 1604 ( 3.8) 609 (10.1) 72 (19.2) 32 ( 50.8)
ventilat (%) 0 23375 (55.2) 1183 (19.6) 68 (18.1) 4 ( 6.3) <0.001
1 18963 (44.8) 4864 (80.4) 307 (81.9) 59 ( 93.7)
mort28 (%) 0 36554 (86.3) 5311 (87.8) 330 (88.0) 49 ( 77.8) 0.002
1 5784 (13.7) 736 (12.2) 45 (12.0) 14 ( 22.2)
mort90 (%) 0 34074 (80.5) 4997 (82.6) 285 (76.0) 35 ( 55.6) <0.001
1 8264 (19.5) 1050 (17.4) 90 (24.0) 28 ( 44.4)
mort365 (%) 0 30331 (71.6) 4627 (76.5) 247 (65.9) 25 ( 39.7) <0.001
1 12007 (28.4) 1420 (23.5) 128 (34.1) 38 ( 60.3)
survival_day (mean (sd)) 497.37 (730.56) 629.44 (862.48) 611.70 (852.84) 244.67 (462.24) <0.001
icu_los (mean (sd)) 3.71 (5.08) 6.05 (8.62) 11.77 (15.73) 27.20 (24.24) <0.001
hosp_los (mean (sd)) 9.37 (9.70) 13.84 (14.49) 22.69 (22.43) 39.78 (34.60) <0.001
# prop.table(table(main2$rrh, main2$mort28, dnn=c("RRH", "28-day Mortality")), 1)
# prop.table(table(main2$rrh, main2$mort90, dnn=c("RRH", "90-day Mortality")), 1)
# prop.table(table(main2$rrh, main2$mort365, dnn=c("RRH", "365-day Mortality")), 1)

CreateTableOne(vars=c("mort28", "mort90", "mort365"),
               strata="rrh_cat",
               data=main2, test=FALSE)
##                  Stratified by rrh_cat
##                    0            [ 1, 3)      [ 3, 5)     [ 5,11]   
##   n               42338         6047         375         63        
##   mort28 = 1 (%)   5784 (13.7)   736 (12.2)   45 (12.0)  14 (22.2) 
##   mort90 = 1 (%)   8264 (19.5)  1050 (17.4)   90 (24.0)  28 (44.4) 
##   mort365 = 1 (%) 12007 (28.4)  1420 (23.5)  128 (34.1)  38 (60.3)
# par(mfrow=c(1, 3))
# plot(main2$rrh_cat, main2$mort28, xlab="RRH", ylab="28-day Mortality")
# plot(main2$rrh_cat, main2$mort90, xlab="RRH", ylab="90-day Mortality")
# plot(main2$rrh_cat, main2$mort365, xlab="RRH", ylab="365-day Mortality")
# 
# p1 <- plot_prop_by_level(main2, "rrh_cat", "mort28", factor.var2="dm")
# p2 <- plot_prop_by_level(main2, "rrh_cat", "mort90", factor.var2="dm")
# p3 <- plot_prop_by_level(main2, "rrh_cat", "mort365", factor.var2="dm")
# grid.arrange(p1, p2, p3, ncol=3)

o1 <- plot_OR_by_level(main2, "rrh_cat", "mort28", include.ref.group.effect=TRUE)
o2 <- plot_OR_by_level(main2, "rrh_cat", "mort90", include.ref.group.effect=TRUE)
o3 <- plot_OR_by_level(main2, "rrh_cat", "mort365", include.ref.group.effect=TRUE)
grid.arrange(o1, o2, o3, ncol=3)

summary(fit.glm <- glm(mort28 ~ rrh, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort28 ~ rrh, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.5420  -0.5420  -0.5420  -0.5114   2.0489  
## 
## Coefficients:
##             Estimate Std. Error  z value Pr(>|z|)    
## (Intercept) -1.84370    0.01415 -130.288  < 2e-16 ***
## rrh1        -0.12443    0.04042   -3.078  0.00208 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 38602  on 48822  degrees of freedom
## Residual deviance: 38592  on 48821  degrees of freedom
## AIC: 38596
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                 2.5 %    97.5 %
## (Intercept) 0.1538889 0.1626666
## rrh1        0.8152508 0.9552405
sjp.glm(fit.glm)
## Waiting for profiling to be done...

summary(fit.glm <- glm(mort90 ~ rrh, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort90 ~ rrh, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.6590  -0.6590  -0.6590  -0.6302   1.8516  
## 
## Coefficients:
##             Estimate Std. Error  z value Pr(>|z|)    
## (Intercept) -1.41663    0.01226 -115.531  < 2e-16 ***
## rrh1        -0.09899    0.03456   -2.864  0.00418 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 47926  on 48822  degrees of freedom
## Residual deviance: 47918  on 48821  degrees of freedom
## AIC: 47922
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                 2.5 %    97.5 %
## (Intercept) 0.2367577 0.2484158
## rrh1        0.8461104 0.9688845
sjp.glm(fit.glm)
## Waiting for profiling to be done...

summary(fit.glm <- glm(mort365 ~ rrh, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort365 ~ rrh, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.8167  -0.8167  -0.8167   1.5876   1.6783  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.92668    0.01078 -85.946  < 2e-16 ***
## rrh1        -0.20114    0.03084  -6.523 6.91e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 57753  on 48822  degrees of freedom
## Residual deviance: 57709  on 48821  degrees of freedom
## AIC: 57713
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                 2.5 %    97.5 %
## (Intercept) 0.3875751 0.4043073
## rrh1        0.7696542 0.8685507
#sjp.glm(fit.glm)


summary(fit.glm <- glm(mort28 ~ rrh_cat, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort28 ~ rrh_cat, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.7090  -0.5420  -0.5420  -0.5095   2.0593  
## 
## Coefficients:
##                Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)    -1.84370    0.01415 -130.288  < 2e-16 ***
## rrh_cat[ 1, 3) -0.13261    0.04180   -3.173  0.00151 ** 
## rrh_cat[ 3, 5) -0.14874    0.15954   -0.932  0.35119    
## rrh_cat[ 5,11]  0.59093    0.30338    1.948  0.05143 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 38602  on 48822  degrees of freedom
## Residual deviance: 38587  on 48819  degrees of freedom
## AIC: 38595
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                    2.5 %    97.5 %
## (Intercept)    0.1538889 0.1626666
## rrh_cat[ 1, 3) 0.8063788 0.9499684
## rrh_cat[ 3, 5) 0.6222577 1.1647193
## rrh_cat[ 5,11] 0.9594442 3.1819085
sjp.glm(fit.glm)
## Waiting for profiling to be done...

summary(fit.glm <- glm(mort90 ~ rrh_cat, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort90 ~ rrh_cat, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.0842  -0.6590  -0.6590  -0.6176   1.8712  
## 
## Coefficients:
##                Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)    -1.41663    0.01226 -115.531  < 2e-16 ***
## rrh_cat[ 1, 3) -0.14342    0.03610   -3.973 7.08e-05 ***
## rrh_cat[ 3, 5)  0.26395    0.12153    2.172   0.0299 *  
## rrh_cat[ 5,11]  1.19348    0.25384    4.702 2.58e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 47926  on 48822  degrees of freedom
## Residual deviance: 47884  on 48819  degrees of freedom
## AIC: 47892
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                    2.5 %    97.5 %
## (Intercept)    0.2367577 0.2484158
## rrh_cat[ 1, 3) 0.8068694 0.9295173
## rrh_cat[ 3, 5) 1.0208289 1.6447957
## rrh_cat[ 5,11] 1.9913018 5.4139628
sjp.glm(fit.glm)
## Waiting for profiling to be done...

summary(fit.glm <- glm(mort365 ~ rrh_cat, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort365 ~ rrh_cat, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3596  -0.8167  -0.8167   1.5876   1.7023  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -0.92668    0.01078 -85.946  < 2e-16 ***
## rrh_cat[ 1, 3) -0.25457    0.03220  -7.907 2.64e-15 ***
## rrh_cat[ 3, 5)  0.26932    0.10944   2.461   0.0139 *  
## rrh_cat[ 5,11]  1.34539    0.25774   5.220 1.79e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 57753  on 48822  degrees of freedom
## Residual deviance: 57652  on 48819  degrees of freedom
## AIC: 57660
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                    2.5 %    97.5 %
## (Intercept)    0.3875751 0.4043073
## rrh_cat[ 1, 3) 0.7276365 0.8255234
## rrh_cat[ 3, 5) 1.0536091 1.6187237
## rrh_cat[ 5,11] 2.3310295 6.4372281
sjp.glm(fit.glm)
## Waiting for profiling to be done...

Bivariate analysis

RH and DM

  1. Compute the OR of 28/90/365-day mortality between different RH groups AND diabetes(-/+).

    • Conclusion: The RH+/DM- group has increased OR, RH+/DM+ group has decreased OR. Can see this phenomenon also in categorized RH.
    • (in discussion) RH may be ignored in non-DM patients can cause damage
CreateTableOne(vars=c("mort28", "mort90", "mort365"),
               strata=c("dm", "rh_cat"),
               data=main2, test=TRUE)
##                  Stratified by dm:rh_cat
##                   0: 0          1: 0          0:[ 1, 3)    1:[ 1, 3)   
##   n               25926         13451         4442         3990        
##   mort28 = 1 (%)   3385 (13.1)   1872 (13.9)   673 (15.2)   516 (12.9) 
##   mort90 = 1 (%)   4738 (18.3)   2786 (20.7)   906 (20.4)   742 (18.6) 
##   mort365 = 1 (%)  6797 (26.2)   4226 (31.4)  1171 (26.4)  1041 (26.1) 
##                  Stratified by dm:rh_cat
##                   0:[ 3, 5)   1:[ 3, 5)   0:[ 5,19]  1:[ 5,19]   p     
##   n               320         508         67         119               
##   mort28 = 1 (%)   50 (15.6)   50 ( 9.8)  15 (22.4)   18 (15.1)  <0.001
##   mort90 = 1 (%)   94 (29.4)   93 (18.3)  33 (49.3)   40 (33.6)  <0.001
##   mort365 = 1 (%) 122 (38.1)  135 (26.6)  40 (59.7)   61 (51.3)  <0.001
##                  Stratified by dm:rh_cat
##                   test
##   n                   
##   mort28 = 1 (%)      
##   mort90 = 1 (%)      
##   mort365 = 1 (%)
main2$rh <- as.factor(main2$rh)

o1 <- plot_OR_by_level(main2, "dm", "mort28", factor.var2="rh", include.ref.group.effect=TRUE)
o2 <- plot_OR_by_level(main2, "dm", "mort90", factor.var2="rh", include.ref.group.effect=TRUE)
o3 <- plot_OR_by_level(main2, "dm", "mort365", factor.var2="rh", include.ref.group.effect=TRUE)
grid.arrange(o1, o2, o3, ncol=3)

o1 <- plot_OR_by_level(main2, "dm", "mort28", factor.var2="rh_cat", include.ref.group.effect=TRUE)
o2 <- plot_OR_by_level(main2, "dm", "mort90", factor.var2="rh_cat", include.ref.group.effect=TRUE)
o3 <- plot_OR_by_level(main2, "dm", "mort365", factor.var2="rh_cat", include.ref.group.effect=TRUE)
grid.arrange(o1, o2, o3, ncol=3)

summary(fit.glm <- glm(mort28 ~ dm + rh, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort28 ~ dm + rh, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.5505  -0.5371  -0.5344  -0.5344   2.0084  
## 
## Coefficients:
##             Estimate Std. Error  z value Pr(>|z|)    
## (Intercept) -1.87400    0.01759 -106.531   <2e-16 ***
## dm1          0.01074    0.02762    0.389    0.697    
## rh1          0.05306    0.03340    1.589    0.112    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 38602  on 48822  degrees of freedom
## Residual deviance: 38599  on 48820  degrees of freedom
## AIC: 38605
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                 2.5 %    97.5 %
## (Intercept) 0.1482856 0.1588719
## dm1         0.9574531 1.0669250
## rh1         0.9873757 1.1255070
sjp.glm(fit.glm)
## Waiting for profiling to be done...

summary(fit.glm <- glm(mort90 ~ dm + rh, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort90 ~ dm + rh, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.6855  -0.6686  -0.6421  -0.6421   1.8333  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.47429    0.01533 -96.151  < 2e-16 ***
## dm1          0.08986    0.02375   3.783 0.000155 ***
## rh1          0.05587    0.02887   1.935 0.053009 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 47926  on 48822  degrees of freedom
## Residual deviance: 47906  on 48820  degrees of freedom
## AIC: 47912
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                 2.5 %   97.5 %
## (Intercept) 0.2221435 0.235905
## dm1         1.0442027 1.146109
## rh1         0.9990864 1.118821
sjp.glm(fit.glm)
## Waiting for profiling to be done...

summary(fit.glm <- glm(mort365 ~ dm + rh, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort365 ~ dm + rh, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.8548  -0.8306  -0.7872   1.5389   1.6572  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.01293    0.01352 -74.902  < 2e-16 ***
## dm1          0.19409    0.02088   9.295  < 2e-16 ***
## rh1         -0.06821    0.02592  -2.632  0.00849 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 57753  on 48822  degrees of freedom
## Residual deviance: 57664  on 48820  degrees of freedom
## AIC: 57670
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                 2.5 %    97.5 %
## (Intercept) 0.3536360 0.3728888
## dm1         1.1654888 1.2649051
## rh1         0.8876947 0.9826150
sjp.glm(fit.glm)
## Waiting for profiling to be done...

summary(fit.glm <- glm(mort28 ~ dm + rh_cat, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort28 ~ dm + rh_cat, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.6262  -0.5373  -0.5344  -0.5344   2.0592  
## 
## Coefficients:
##               Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)   -1.87429    0.01760 -106.508   <2e-16 ***
## dm1            0.01158    0.02764    0.419   0.6754    
## rh_cat[ 1, 3)  0.06188    0.03481    1.777   0.0755 .  
## rh_cat[ 3, 5) -0.11796    0.10793   -1.093   0.2745    
## rh_cat[ 5,19]  0.33294    0.19268    1.728   0.0840 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 38602  on 48822  degrees of freedom
## Residual deviance: 38594  on 48818  degrees of freedom
## AIC: 38604
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                   2.5 %    97.5 %
## (Intercept)   0.1482412 0.1588282
## dm1           0.9582095 1.0678646
## rh_cat[ 1, 3) 0.9933160 1.1385630
## rh_cat[ 3, 5) 0.7152436 1.0924396
## rh_cat[ 5,19] 0.9407063 2.0071034
sjp.glm(fit.glm)
## Waiting for profiling to be done...

summary(fit.glm <- glm(mort90 ~ dm + rh_cat, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort90 ~ dm + rh_cat, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.0104  -0.6677  -0.6426  -0.6426   1.8325  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -1.47265    0.01533 -96.045  < 2e-16 ***
## dm1            0.08522    0.02380   3.581 0.000342 ***
## rh_cat[ 1, 3)  0.01677    0.03048   0.550 0.582182    
## rh_cat[ 3, 5)  0.18797    0.08435   2.228 0.025858 *  
## rh_cat[ 5,19]  0.98102    0.15090   6.501 7.96e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 47926  on 48822  degrees of freedom
## Residual deviance: 47867  on 48818  degrees of freedom
## AIC: 47877
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                   2.5 %    97.5 %
## (Intercept)   0.2225084 0.2362923
## dm1           1.0392815 1.1409000
## rh_cat[ 1, 3) 0.9577248 1.0792562
## rh_cat[ 3, 5) 1.0203872 1.4205445
## rh_cat[ 5,19] 1.9771240 3.5758771
sjp.glm(fit.glm)
## Waiting for profiling to be done...

summary(fit.glm <- glm(mort365 ~ dm + rh_cat, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort365 ~ dm + rh_cat, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.2813  -0.8134  -0.7877   1.5404   1.6772  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -1.01114    0.01353 -74.760  < 2e-16 ***
## dm1            0.18913    0.02092   9.041  < 2e-16 ***
## rh_cat[ 1, 3) -0.11434    0.02735  -4.180 2.91e-05 ***
## rh_cat[ 3, 5)  0.09517    0.07623   1.248    0.212    
## rh_cat[ 5,19]  1.06293    0.14788   7.188 6.60e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 57753  on 48822  degrees of freedom
## Residual deviance: 57598  on 48818  degrees of freedom
## AIC: 57608
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                   2.5 %    97.5 %
## (Intercept)   0.3542670 0.3735566
## dm1           1.1596372 1.2587446
## rh_cat[ 1, 3) 0.8452661 0.9409359
## rh_cat[ 3, 5) 0.9458791 1.2754553
## rh_cat[ 5,19] 2.1679138 3.8746577
sjp.glm(fit.glm)
## Waiting for profiling to be done...

RRH and DM

  1. Compute the OR of 28/90/365-day mortality between different RRH groups AND diabetes(-/+).

    • Conclusion: The RRH+ group has decreased OR in most cases (except 90-day mortality in non-DM).
    • (in discussion) RRH may increase the attention of the clinicians and maybe they will do something to deal with this rapid glucose level change???
CreateTableOne(vars=c("mort28", "mort90", "mort365"),
               strata=c("dm", "rrh_cat"),
               data=main2, test=FALSE)
##                  Stratified by dm:rrh_cat
##                   0: 0          1: 0          0:[ 1, 3)    1:[ 1, 3)   
##   n               27356         14982         3216         2831        
##   mort28 = 1 (%)   3674 (13.4)   2110 (14.1)   418 (13.0)   318 (11.2) 
##   mort90 = 1 (%)   5125 (18.7)   3139 (21.0)   583 (18.1)   467 (16.5) 
##   mort365 = 1 (%)  7287 (26.6)   4720 (31.5)   764 (23.8)   656 (23.2) 
##                  Stratified by dm:rrh_cat
##                   0:[ 3, 5)   1:[ 3, 5)   0:[ 5,11]  1:[ 5,11] 
##   n               152         223         31         32        
##   mort28 = 1 (%)   23 (15.1)   22 ( 9.9)   8 (25.8)   6 (18.8) 
##   mort90 = 1 (%)   49 (32.2)   41 (18.4)  14 (45.2)  14 (43.8) 
##   mort365 = 1 (%)  60 (39.5)   68 (30.5)  19 (61.3)  19 (59.4)
main2$rrh <- as.factor(main2$rrh)
o1 <- plot_OR_by_level(main2, "dm", "mort28", factor.var2="rrh", include.ref.group.effect=TRUE)
o2 <- plot_OR_by_level(main2, "dm", "mort90", factor.var2="rrh", include.ref.group.effect=TRUE)
o3 <- plot_OR_by_level(main2, "dm", "mort365", factor.var2="rrh", include.ref.group.effect=TRUE)
grid.arrange(o1, o2, o3, ncol=3)

o1 <- plot_OR_by_level(main2, "dm", "mort28", factor.var2="rrh_cat", include.ref.group.effect=TRUE)
o2 <- plot_OR_by_level(main2, "dm", "mort90", factor.var2="rrh_cat", include.ref.group.effect=TRUE)
o3 <- plot_OR_by_level(main2, "dm", "mort365", factor.var2="rrh_cat", include.ref.group.effect=TRUE)
grid.arrange(o1, o2, o3, ncol=3)

summary(fit.glm <- glm(mort28 ~ dm + rrh, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort28 ~ dm + rrh, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.5458  -0.5458  -0.5399  -0.5144   2.0537  
## 
## Coefficients:
##             Estimate Std. Error  z value Pr(>|z|)    
## (Intercept) -1.85204    0.01724 -107.447   <2e-16 ***
## dm1          0.02345    0.02751    0.852   0.3940    
## rrh1        -0.12730    0.04056   -3.138   0.0017 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 38602  on 48822  degrees of freedom
## Residual deviance: 38591  on 48820  degrees of freedom
## AIC: 38597
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                 2.5 %    97.5 %
## (Intercept) 0.1516841 0.1622874
## dm1         0.9699056 1.0803425
## rrh1        0.8126961 0.9527719
sjp.glm(fit.glm)
## Waiting for profiling to be done...

summary(fit.glm <- glm(mort90 ~ dm + rrh, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort90 ~ dm + rrh, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.6785  -0.6482  -0.6482  -0.6454   1.8734  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.45344    0.01504 -96.619  < 2e-16 ***
## dm1          0.10199    0.02366   4.310 1.63e-05 ***
## rrh1        -0.11154    0.03469  -3.215   0.0013 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 47926  on 48822  degrees of freedom
## Residual deviance: 47899  on 48820  degrees of freedom
## AIC: 47905
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                 2.5 %    97.5 %
## (Intercept) 0.2269527 0.2407383
## dm1         1.0571362 1.1598941
## rrh1        0.8353468 0.9570504
sjp.glm(fit.glm)
## Waiting for profiling to be done...

summary(fit.glm <- glm(mort365 ~ dm + rrh, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort365 ~ dm + rrh, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.8619  -0.7916  -0.7916   1.5299   1.7224  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.99970    0.01330 -75.162  < 2e-16 ***
## dm1          0.20077    0.02082   9.645  < 2e-16 ***
## rrh1        -0.22622    0.03098  -7.302 2.83e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 57753  on 48822  degrees of freedom
## Residual deviance: 57616  on 48820  degrees of freedom
## AIC: 57622
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                 2.5 %    97.5 %
## (Intercept) 0.3585030 0.3776905
## dm1         1.1734476 1.2732153
## rrh1        0.7503759 0.8472720
sjp.glm(fit.glm)
## Waiting for profiling to be done...

summary(fit.glm <- glm(mort28 ~ dm + rrh_cat, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort28 ~ dm + rrh_cat, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.7126  -0.5458  -0.5399  -0.5124   2.0652  
## 
## Coefficients:
##                Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)    -1.85201    0.01724 -107.436  < 2e-16 ***
## dm1             0.02337    0.02752    0.849  0.39578    
## rrh_cat[ 1, 3) -0.13529    0.04192   -3.227  0.00125 ** 
## rrh_cat[ 3, 5) -0.15437    0.15968   -0.967  0.33367    
## rrh_cat[ 5,11]  0.58734    0.30341    1.936  0.05289 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 38602  on 48822  degrees of freedom
## Residual deviance: 38586  on 48818  degrees of freedom
## AIC: 38596
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                    2.5 %    97.5 %
## (Intercept)    0.1516881 0.1622926
## dm1            0.9698101 1.0802705
## rrh_cat[ 1, 3) 0.8040375 0.9476541
## rrh_cat[ 3, 5) 0.6186078 1.1585182
## rrh_cat[ 5,11] 0.9559488 3.1707204
sjp.glm(fit.glm)
## Waiting for profiling to be done...

summary(fit.glm <- glm(mort90 ~ dm + rrh_cat, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort90 ~ dm + rrh_cat, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.1045  -0.6484  -0.6484  -0.6325   1.8923  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -1.45277    0.01504 -96.572  < 2e-16 ***
## dm1             0.10017    0.02368   4.229 2.34e-05 ***
## rrh_cat[ 1, 3) -0.15499    0.03621  -4.281 1.86e-05 ***
## rrh_cat[ 3, 5)  0.23990    0.12169   1.971   0.0487 *  
## rrh_cat[ 5,11]  1.17861    0.25394   4.641 3.46e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 47926  on 48822  degrees of freedom
## Residual deviance: 47866  on 48818  degrees of freedom
## AIC: 47876
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                    2.5 %    97.5 %
## (Intercept)    0.2271050 0.2409003
## dm1            1.0551684 1.1578213
## rrh_cat[ 1, 3) 0.7974180 0.9190325
## rrh_cat[ 3, 5) 0.9962721 1.6062249
## rrh_cat[ 5,11] 1.9615259 5.3350886
sjp.glm(fit.glm)
## Waiting for profiling to be done...

summary(fit.glm <- glm(mort365 ~ dm + rrh_cat, data=main2, family="binomial"))
## 
## Call:
## glm(formula = mort365 ~ dm + rrh_cat, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.4036  -0.7919  -0.7919   1.5305   1.7452  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -0.99894    0.01330 -75.095  < 2e-16 ***
## dm1             0.19874    0.02084   9.538  < 2e-16 ***
## rrh_cat[ 1, 3) -0.27796    0.03232  -8.599  < 2e-16 ***
## rrh_cat[ 3, 5)  0.22187    0.10967   2.023   0.0431 *  
## rrh_cat[ 5,11]  1.31772    0.25806   5.106 3.29e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 57753  on 48822  degrees of freedom
## Residual deviance: 57561  on 48818  degrees of freedom
## AIC: 57571
## 
## Number of Fisher Scoring iterations: 4
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                    2.5 %    97.5 %
## (Intercept)    0.3587724 0.3779770
## dm1            1.1710320 1.2706937
## rrh_cat[ 1, 3) 0.7106387 0.8066419
## rrh_cat[ 3, 5) 1.0043393 1.5443958
## rrh_cat[ 5,11] 2.2659738 6.2652578
sjp.glm(fit.glm)
## Waiting for profiling to be done...

Multivariate analysis

insNo <- main2[main2$ins==0, ]
insIv <- main2[main2$ins==1 & main2$ins_iv==1 & main2$ins_sc==0, ]
insSc <- main2[main2$ins==1 & main2$ins_iv==0 & main2$ins_sc==1, ]
insIS <- main2[main2$ins==1 & main2$ins_iv==1 & main2$ins_sc==1, ]

RH + diabetes

No insulin

# no insulin group
o1 <- plot_OR_by_level(insNo, "dm", "mort28", factor.var2="rh", include.ref.group.effect=TRUE)
o2 <- plot_OR_by_level(insNo, "dm", "mort90", factor.var2="rh", include.ref.group.effect=TRUE)
o3 <- plot_OR_by_level(insNo, "dm", "mort365", factor.var2="rh", include.ref.group.effect=TRUE)
grid.arrange(o1, o2, o3, ncol=3)

o1 <- plot_OR_by_level(insNo, "dm", "mort28", factor.var2="rh_cat", include.ref.group.effect=TRUE)
o2 <- plot_OR_by_level(insNo, "dm", "mort90", factor.var2="rh_cat", include.ref.group.effect=TRUE)
o3 <- plot_OR_by_level(insNo, "dm", "mort365", factor.var2="rh_cat", include.ref.group.effect=TRUE)
grid.arrange(o1, o2, o3, ncol=3)
## Warning: Removed 1 rows containing missing values (geom_errorbar).

## Warning: Removed 1 rows containing missing values (geom_errorbar).

## Warning: Removed 1 rows containing missing values (geom_errorbar).

IV insulin

# IV insulin group
o1 <- plot_OR_by_level(insIv, "dm", "mort28", factor.var2="rh", include.ref.group.effect=TRUE)
o2 <- plot_OR_by_level(insIv, "dm", "mort90", factor.var2="rh", include.ref.group.effect=TRUE)
o3 <- plot_OR_by_level(insIv, "dm", "mort365", factor.var2="rh", include.ref.group.effect=TRUE)
grid.arrange(o1, o2, o3, ncol=3)

o1 <- plot_OR_by_level(insIv, "dm", "mort28", factor.var2="rh_cat", include.ref.group.effect=TRUE)
o2 <- plot_OR_by_level(insIv, "dm", "mort90", factor.var2="rh_cat", include.ref.group.effect=TRUE)
o3 <- plot_OR_by_level(insIv, "dm", "mort365", factor.var2="rh_cat", include.ref.group.effect=TRUE)
grid.arrange(o1, o2, o3, ncol=3)

Subcutaneous insulin

# subcutaneous insulin group
o1 <- plot_OR_by_level(insSc, "dm", "mort28", factor.var2="rh", include.ref.group.effect=TRUE)
o2 <- plot_OR_by_level(insSc, "dm", "mort90", factor.var2="rh", include.ref.group.effect=TRUE)
o3 <- plot_OR_by_level(insSc, "dm", "mort365", factor.var2="rh", include.ref.group.effect=TRUE)
grid.arrange(o1, o2, o3, ncol=3)

o1 <- plot_OR_by_level(insSc, "dm", "mort28", factor.var2="rh_cat", include.ref.group.effect=TRUE)
o2 <- plot_OR_by_level(insSc, "dm", "mort90", factor.var2="rh_cat", include.ref.group.effect=TRUE)
o3 <- plot_OR_by_level(insSc, "dm", "mort365", factor.var2="rh_cat", include.ref.group.effect=TRUE)
grid.arrange(o1, o2, o3, ncol=3)
## Warning: Removed 3 rows containing missing values (geom_errorbar).
## Warning: Removed 2 rows containing missing values (geom_errorbar).

## Warning: Removed 2 rows containing missing values (geom_errorbar).

IV + subcutaneous insulin

# IV+subcutaneous insulin group
o1 <- plot_OR_by_level(insIS, "dm", "mort28", factor.var2="rh", include.ref.group.effect=TRUE)
o2 <- plot_OR_by_level(insIS, "dm", "mort90", factor.var2="rh", include.ref.group.effect=TRUE)
o3 <- plot_OR_by_level(insIS, "dm", "mort365", factor.var2="rh", include.ref.group.effect=TRUE)
grid.arrange(o1, o2, o3, ncol=3)

o1 <- plot_OR_by_level(insIS, "dm", "mort28", factor.var2="rh_cat", include.ref.group.effect=TRUE)
o2 <- plot_OR_by_level(insIS, "dm", "mort90", factor.var2="rh_cat", include.ref.group.effect=TRUE)
o3 <- plot_OR_by_level(insIS, "dm", "mort365", factor.var2="rh_cat", include.ref.group.effect=TRUE)
grid.arrange(o1, o2, o3, ncol=3)

RRH

No insulin

# no insulin group
o1 <- plot_OR_by_level(insNo, "dm", "mort28", factor.var2="rrh", include.ref.group.effect=TRUE)
o2 <- plot_OR_by_level(insNo, "dm", "mort90", factor.var2="rrh", include.ref.group.effect=TRUE)
o3 <- plot_OR_by_level(insNo, "dm", "mort365", factor.var2="rrh", include.ref.group.effect=TRUE)
grid.arrange(o1, o2, o3, ncol=3)

o1 <- plot_OR_by_level(insNo, "dm", "mort28", factor.var2="rrh_cat", include.ref.group.effect=TRUE)
o2 <- plot_OR_by_level(insNo, "dm", "mort90", factor.var2="rrh_cat", include.ref.group.effect=TRUE)
o3 <- plot_OR_by_level(insNo, "dm", "mort365", factor.var2="rrh_cat", include.ref.group.effect=TRUE)
grid.arrange(o1, o2, o3, ncol=3)
## Warning: Removed 1 rows containing missing values (geom_errorbar).

## Warning: Removed 1 rows containing missing values (geom_errorbar).

## Warning: Removed 1 rows containing missing values (geom_errorbar).

IV insulin

# IV insulin group
o1 <- plot_OR_by_level(insIv, "dm", "mort28", factor.var2="rrh", include.ref.group.effect=TRUE)
o2 <- plot_OR_by_level(insIv, "dm", "mort90", factor.var2="rrh", include.ref.group.effect=TRUE)
o3 <- plot_OR_by_level(insIv, "dm", "mort365", factor.var2="rrh", include.ref.group.effect=TRUE)
grid.arrange(o1, o2, o3, ncol=3)

o1 <- plot_OR_by_level(insIv, "dm", "mort28", factor.var2="rrh_cat", include.ref.group.effect=TRUE)
o2 <- plot_OR_by_level(insIv, "dm", "mort90", factor.var2="rrh_cat", include.ref.group.effect=TRUE)
o3 <- plot_OR_by_level(insIv, "dm", "mort365", factor.var2="rrh_cat", include.ref.group.effect=TRUE)
grid.arrange(o1, o2, o3, ncol=3)

Subcutaneous insulin

# subcutaneous insulin group
o1 <- plot_OR_by_level(insSc, "dm", "mort28", factor.var2="rrh", include.ref.group.effect=TRUE)
o2 <- plot_OR_by_level(insSc, "dm", "mort90", factor.var2="rrh", include.ref.group.effect=TRUE)
o3 <- plot_OR_by_level(insSc, "dm", "mort365", factor.var2="rrh", include.ref.group.effect=TRUE)
grid.arrange(o1, o2, o3, ncol=3)

o1 <- plot_OR_by_level(insSc, "dm", "mort28", factor.var2="rrh_cat", include.ref.group.effect=TRUE)
o2 <- plot_OR_by_level(insSc, "dm", "mort90", factor.var2="rrh_cat", include.ref.group.effect=TRUE)
o3 <- plot_OR_by_level(insSc, "dm", "mort365", factor.var2="rrh_cat", include.ref.group.effect=TRUE)
grid.arrange(o1, o2, o3, ncol=3)
## Warning: Removed 3 rows containing missing values (geom_errorbar).

## Warning: Removed 3 rows containing missing values (geom_errorbar).

## Warning: Removed 3 rows containing missing values (geom_errorbar).

IV + subcutaneous insulin

# IV+subcutaneous insulin group
o1 <- plot_OR_by_level(insIS, "dm", "mort28", factor.var2="rrh", include.ref.group.effect=TRUE)
o2 <- plot_OR_by_level(insIS, "dm", "mort90", factor.var2="rrh", include.ref.group.effect=TRUE)
o3 <- plot_OR_by_level(insIS, "dm", "mort365", factor.var2="rrh", include.ref.group.effect=TRUE)
grid.arrange(o1, o2, o3, ncol=3)

o1 <- plot_OR_by_level(insIS, "dm", "mort28", factor.var2="rrh_cat", include.ref.group.effect=TRUE)
o2 <- plot_OR_by_level(insIS, "dm", "mort90", factor.var2="rrh_cat", include.ref.group.effect=TRUE)
o3 <- plot_OR_by_level(insIS, "dm", "mort365", factor.var2="rrh_cat", include.ref.group.effect=TRUE)
grid.arrange(o1, o2, o3, ncol=3)
## Warning: Removed 1 rows containing missing values (geom_errorbar).

END

Propensity score model (skip)

  1. We used propensity score model to adjust the confounding factors, which include age, gender, SOFA score, Elixhauser morbidity index, insulin, TPN, vasopressor and mechanical ventilator use.
  2. The age and SOFA scores were categorized into four groups, and the insulin, TPN, vasopressor and mechanical ventilator were transformed into binary variables.
  3. Stochastic gradient boosting algorithm was used to find the optimal prediction model to balance the RH and non-RH groups.

    • Model 1: RH ~ age+gender+sofa+elixhauser+insulin+TPN+vasopressor+ventilator
    • Model 2: RRH ~ age+gender+sofa+elixhauser+insulin+TPN+vasopressor+ventilator
  4. The propsensity score was stratified to four groups. The admission with higher propensity score indicates the higher probability of having RH.
  5. However, the propensity score groups (for both RH and RRH) seem not balanced after using the model. (why??????)

    • Conclusion: In general, OR is higher in moderate propensity score (both RH and RRH), and higher in non-diabetes -> usually be ignored???

Model 1: RH

Model 2: RRH

We also constructed the propensity score model for RRH.

Logistic regression (the result of this part is quite confusing)

Next, we created the logistic regression model for 28-day mortality based on the propensity score of RH in diabetes and non-diabetes groups.

  1. I fit two logistic regression models:

    • Model 1: diabetes + propensity score (skip)
    • Model 2: diabetes + RH + all other covariates (didn’t use propensity score)
  2. If i use model1 (skip): propensity score is significant to increase the OR, but diabetes is not significant (but decrease the OR) with p-value around 0.06. Non-diabetes and diabetes have similar trend.

  3. If i use all variables, with binary/categorized RH: diabetes(+) increases OR, RH(+) decreases OR -> this one is not agree with the previous observation (inverse)

# fit.lm <- lm(mort28 ~ rh + dm + ins_mean + age + gender + sofa_cat + elixhauser_hospital + tpn + vasop + ventilat, 
#              data=main2)
# summary(fit.lm)
# confint(fit.lm)
# sjp.lm(fit.lm, type="coef")

fit.glm <- glm(mort28 ~ rh + dm + ins_iv + ins_sc + age + gender + sofa_1 + elixhauser_hospital + tpn + vasop + ventilat,
               data=main2, 
               family="binomial")
summary(fit.glm)
## 
## Call:
## glm(formula = mort28 ~ rh + dm + ins_iv + ins_sc + age + gender + 
##     sofa_1 + elixhauser_hospital + tpn + vasop + ventilat, family = "binomial", 
##     data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5229  -0.5301  -0.3710  -0.2340   3.0847  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -4.7992383  0.0800406 -59.960  < 2e-16 ***
## rh1                 -0.3089441  0.0403406  -7.658 1.88e-14 ***
## dm1                  0.1390234  0.0308431   4.507 6.56e-06 ***
## ins_iv1             -0.5451002  0.0348841 -15.626  < 2e-16 ***
## ins_sc1             -0.4185274  0.0338552 -12.362  < 2e-16 ***
## age                  0.0314106  0.0009989  31.445  < 2e-16 ***
## gender1             -0.1237774  0.0294055  -4.209 2.56e-05 ***
## sofa_1               0.1450846  0.0063289  22.924  < 2e-16 ***
## elixhauser_hospital  0.0800924  0.0023626  33.900  < 2e-16 ***
## tpn1                -0.0884337  0.0608283  -1.454    0.146    
## vasop1               1.8005852  0.0520050  34.623  < 2e-16 ***
## ventilat1            0.4175757  0.0380739  10.967  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 38569  on 48813  degrees of freedom
## Residual deviance: 32005  on 48802  degrees of freedom
##   (9 observations deleted due to missingness)
## AIC: 32029
## 
## Number of Fisher Scoring iterations: 5
exp(coef(fit.glm)[1])
## (Intercept) 
## 0.008236018
#confint(fit.glm)
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                           2.5 %      97.5 %
## (Intercept)         0.007033838 0.009626535
## rh1                 0.678186691 0.794387027
## dm1                 1.081663548 1.220684849
## ins_iv1             0.541357949 0.620692275
## ins_sc1             0.615646743 0.703026559
## age                 1.029897881 1.033938935
## gender1             0.834102825 0.936016156
## sofa_1              1.141896690 1.170582342
## elixhauser_hospital 1.078387298 1.088421550
## tpn1                0.811797029 1.030430821
## vasop1              5.467220853 6.703626045
## ventilat1           1.409090932 1.635915984
sjp.glm(fit.glm)
## Waiting for profiling to be done...

fit.glm <- glm(mort28 ~ rh_freq + dm + ins_iv + ins_sc + age + gender + sofa_1 + elixhauser_hospital + tpn + vasop + ventilat,
               data=main2, 
               family="binomial")
summary(fit.glm)
## 
## Call:
## glm(formula = mort28 ~ rh_freq + dm + ins_iv + ins_sc + age + 
##     gender + sofa_1 + elixhauser_hospital + tpn + vasop + ventilat, 
##     family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4701  -0.5300  -0.3712  -0.2341   3.0913  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -4.8098089  0.0801284 -60.026  < 2e-16 ***
## rh_freq             -0.1839124  0.0214360  -8.580  < 2e-16 ***
## dm1                  0.1428346  0.0308590   4.629 3.68e-06 ***
## ins_iv1             -0.5330194  0.0348838 -15.280  < 2e-16 ***
## ins_sc1             -0.4247263  0.0338968 -12.530  < 2e-16 ***
## age                  0.0314572  0.0009999  31.462  < 2e-16 ***
## gender1             -0.1261649  0.0294268  -4.287 1.81e-05 ***
## sofa_1               0.1446291  0.0063336  22.835  < 2e-16 ***
## elixhauser_hospital  0.0801042  0.0023624  33.908  < 2e-16 ***
## tpn1                -0.0488845  0.0612317  -0.798    0.425    
## vasop1               1.8283346  0.0524956  34.828  < 2e-16 ***
## ventilat1            0.4113932  0.0380070  10.824  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 38569  on 48813  degrees of freedom
## Residual deviance: 31983  on 48802  degrees of freedom
##   (9 observations deleted due to missingness)
## AIC: 32007
## 
## Number of Fisher Scoring iterations: 5
exp(coef(fit.glm)[1])
## (Intercept) 
## 0.008149417
#confint(fit.glm)
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                           2.5 %      97.5 %
## (Intercept)         0.006958664 0.009526932
## rh_freq             0.797353310 0.867265773
## dm1                 1.085760425 1.225384746
## ins_iv1             0.547937049 0.628234929
## ins_sc1             0.611791729 0.698738335
## age                 1.029943941 1.033988995
## gender1             0.832078808 0.933822800
## sofa_1              1.141366245 1.170059896
## elixhauser_hospital 1.078400494 1.088434031
## tpn1                0.843891668 1.072864768
## vasop1              5.615797123 6.899066018
## ventilat1           1.400591096 1.625621422
sjp.glm(fit.glm)
## Waiting for profiling to be done...

  1. Remove TPN, vasopressor, ventilator, with binary/categorized RH: the same
fit.glm <- glm(mort28 ~ rh + dm + ins_iv + ins_sc + age + gender + sofa_1 + elixhauser_hospital,
               data=main2, 
               family="binomial")
summary(fit.glm)
## 
## Call:
## glm(formula = mort28 ~ rh + dm + ins_iv + ins_sc + age + gender + 
##     sofa_1 + elixhauser_hospital, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1840  -0.5559  -0.3886  -0.2447   3.0492  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -4.6639762  0.0769785 -60.588  < 2e-16 ***
## rh1                 -0.1487115  0.0380070  -3.913 9.13e-05 ***
## dm1                  0.1445066  0.0300511   4.809 1.52e-06 ***
## ins_iv1             -0.4355772  0.0334276 -13.030  < 2e-16 ***
## ins_sc1             -0.3721726  0.0327801 -11.354  < 2e-16 ***
## age                  0.0295455  0.0009662  30.580  < 2e-16 ***
## gender1             -0.1293699  0.0286964  -4.508 6.54e-06 ***
## sofa_1               0.2232781  0.0051158  43.645  < 2e-16 ***
## elixhauser_hospital  0.0843759  0.0022781  37.038  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 38569  on 48813  degrees of freedom
## Residual deviance: 33327  on 48805  degrees of freedom
##   (9 observations deleted due to missingness)
## AIC: 33345
## 
## Number of Fisher Scoring iterations: 5
exp(coef(fit.glm)[1])
## (Intercept) 
## 0.009428897
#confint(fit.glm)
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                          2.5 %     97.5 %
## (Intercept)         0.00810154 0.01095535
## rh1                 0.79973717 0.92823001
## dm1                 1.08930813 1.22549869
## ins_iv1             0.60576825 0.69058519
## ins_sc1             0.64622889 0.73484369
## age                 1.02804429 1.03194545
## gender1             0.83060634 0.92950299
## sofa_1              1.23771568 1.26278860
## elixhauser_hospital 1.08319655 1.09291331
sjp.glm(fit.glm)
## Waiting for profiling to be done...

fit.glm <- glm(mort28 ~ rh_freq + dm + ins_iv + ins_sc + age + gender + sofa_1 + elixhauser_hospital,
               data=main2, 
               family="binomial")
summary(fit.glm)
## 
## Call:
## glm(formula = mort28 ~ rh_freq + dm + ins_iv + ins_sc + age + 
##     gender + sofa_1 + elixhauser_hospital, family = "binomial", 
##     data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1505  -0.5563  -0.3887  -0.2445   3.0506  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -4.6680280  0.0769664 -60.650  < 2e-16 ***
## rh_freq             -0.0787647  0.0199821  -3.942 8.09e-05 ***
## dm1                  0.1454696  0.0300568   4.840 1.30e-06 ***
## ins_iv1             -0.4335937  0.0334680 -12.955  < 2e-16 ***
## ins_sc1             -0.3750104  0.0327968 -11.434  < 2e-16 ***
## age                  0.0295331  0.0009663  30.564  < 2e-16 ***
## gender1             -0.1306100  0.0287045  -4.550 5.36e-06 ***
## sofa_1               0.2229119  0.0050956  43.746  < 2e-16 ***
## elixhauser_hospital  0.0845442  0.0022774  37.123  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 38569  on 48813  degrees of freedom
## Residual deviance: 33327  on 48805  degrees of freedom
##   (9 observations deleted due to missingness)
## AIC: 33345
## 
## Number of Fisher Scoring iterations: 5
exp(coef(fit.glm)[1])
## (Intercept) 
##  0.00939077
#confint(fit.glm)
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                           2.5 %     97.5 %
## (Intercept)         0.008068967 0.01091079
## rh_freq             0.888302618 0.96069656
## dm1                 1.090345258 1.22669285
## ins_iv1             0.606919836 0.69200758
## ins_sc1             0.644376415 0.73278531
## age                 1.028031275 1.03193285
## gender1             0.829563470 0.92836572
## sofa_1              1.237311836 1.26227626
## elixhauser_hospital 1.083380243 1.09309580
sjp.glm(fit.glm)
## Waiting for profiling to be done...

  1. RRH the same
fit.glm <- glm(mort28 ~ rrh + dm + ins_iv + ins_sc + age + gender + sofa_1 + elixhauser_hospital + tpn + vasop + ventilat,
               data=main2, 
               family="binomial")
summary(fit.glm)
## 
## Call:
## glm(formula = mort28 ~ rrh + dm + ins_iv + ins_sc + age + gender + 
##     sofa_1 + elixhauser_hospital + tpn + vasop + ventilat, family = "binomial", 
##     data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6367  -0.5299  -0.3698  -0.2338   3.0919  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -4.7974382  0.0800740 -59.913  < 2e-16 ***
## rrh1                -0.4576246  0.0480232  -9.529  < 2e-16 ***
## dm1                  0.1334795  0.0308276   4.330 1.49e-05 ***
## ins_iv1             -0.5367235  0.0346938 -15.470  < 2e-16 ***
## ins_sc1             -0.4207960  0.0338836 -12.419  < 2e-16 ***
## age                  0.0313867  0.0009995  31.403  < 2e-16 ***
## gender1             -0.1254764  0.0294261  -4.264 2.01e-05 ***
## sofa_1               0.1451927  0.0063238  22.960  < 2e-16 ***
## elixhauser_hospital  0.0793630  0.0023642  33.568  < 2e-16 ***
## tpn1                -0.0802169  0.0608115  -1.319    0.187    
## vasop1               1.8039388  0.0520585  34.652  < 2e-16 ***
## ventilat1            0.4185573  0.0379994  11.015  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 38569  on 48813  degrees of freedom
## Residual deviance: 31969  on 48802  degrees of freedom
##   (9 observations deleted due to missingness)
## AIC: 31993
## 
## Number of Fisher Scoring iterations: 5
exp(coef(fit.glm)[1])
## (Intercept) 
## 0.008250857
#confint(fit.glm)
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                           2.5 %      97.5 %
## (Intercept)         0.007046039 0.009644498
## rrh1                0.575605981 0.694850880
## dm1                 1.075715254 1.213898462
## ins_iv1             0.546115524 0.625680198
## ins_sc1             0.614217215 0.701472258
## age                 1.029872048 1.033915309
## gender1             0.832652833 0.934464799
## sofa_1              1.142031387 1.170697026
## elixhauser_hospital 1.077597495 1.087631452
## tpn1                0.818521307 1.038897992
## vasop1              5.485034940 6.726882382
## ventilat1           1.410680954 1.637283514
sjp.glm(fit.glm)
## Waiting for profiling to be done...

fit.glm <- glm(mort28 ~ rrh_freq + dm + ins_iv + ins_sc + age + gender + sofa_1 + elixhauser_hospital + tpn + vasop + ventilat,
               data=main2, 
               family="binomial")
summary(fit.glm)
## 
## Call:
## glm(formula = mort28 ~ rrh_freq + dm + ins_iv + ins_sc + age + 
##     gender + sofa_1 + elixhauser_hospital + tpn + vasop + ventilat, 
##     family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6323  -0.5297  -0.3704  -0.2340   3.0997  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -4.8030639  0.0800974 -59.965  < 2e-16 ***
## rrh_freq            -0.2697498  0.0306965  -8.788  < 2e-16 ***
## dm1                  0.1340147  0.0308267   4.347 1.38e-05 ***
## ins_iv1             -0.5406248  0.0346828 -15.588  < 2e-16 ***
## ins_sc1             -0.4246465  0.0338989 -12.527  < 2e-16 ***
## age                  0.0314178  0.0009997  31.427  < 2e-16 ***
## gender1             -0.1263220  0.0294275  -4.293 1.77e-05 ***
## sofa_1               0.1443560  0.0063256  22.821  < 2e-16 ***
## elixhauser_hospital  0.0795823  0.0023634  33.673  < 2e-16 ***
## tpn1                -0.0676232  0.0609528  -1.109    0.267    
## vasop1               1.8124414  0.0522114  34.714  < 2e-16 ***
## ventilat1            0.4106044  0.0379768  10.812  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 38569  on 48813  degrees of freedom
## Residual deviance: 31979  on 48802  degrees of freedom
##   (9 observations deleted due to missingness)
## AIC: 32003
## 
## Number of Fisher Scoring iterations: 5
exp(coef(fit.glm)[1])
## (Intercept) 
## 0.008204571
#confint(fit.glm)
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                           2.5 %      97.5 %
## (Intercept)         0.007006186 0.009590827
## rrh_freq            0.718450486 0.810344418
## dm1                 1.076293121 1.214546248
## ins_iv1             0.544000171 0.623229652
## ins_sc1             0.611837920 0.698796842
## age                 1.029903631 1.033947882
## gender1             0.831946721 0.933677411
## sofa_1              1.141072287 1.169721952
## elixhauser_hospital 1.077835677 1.087868129
## tpn1                0.828670400 1.052362530
## vasop1              5.530265316 6.786421337
## ventilat1           1.399568852 1.624242673
sjp.glm(fit.glm)
## Waiting for profiling to be done...

fit.glm <- glm(mort28 ~ rrh + dm + ins_iv + ins_sc + age + gender + sofa_1 + elixhauser_hospital,
               data=main2, 
               family="binomial")
summary(fit.glm)
## 
## Call:
## glm(formula = mort28 ~ rrh + dm + ins_iv + ins_sc + age + gender + 
##     sofa_1 + elixhauser_hospital, family = "binomial", data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2528  -0.5554  -0.3880  -0.2449   3.0808  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -4.6587957  0.0769862 -60.515  < 2e-16 ***
## rrh1                -0.2940796  0.0452914  -6.493 8.41e-11 ***
## dm1                  0.1427694  0.0300296   4.754 1.99e-06 ***
## ins_iv1             -0.4193475  0.0332232 -12.622  < 2e-16 ***
## ins_sc1             -0.3745482  0.0328086 -11.416  < 2e-16 ***
## age                  0.0295009  0.0009664  30.525  < 2e-16 ***
## gender1             -0.1315777  0.0287134  -4.582 4.60e-06 ***
## sofa_1               0.2248910  0.0050866  44.212  < 2e-16 ***
## elixhauser_hospital  0.0838595  0.0022796  36.787  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 38569  on 48813  degrees of freedom
## Residual deviance: 33299  on 48805  degrees of freedom
##   (9 observations deleted due to missingness)
## AIC: 33317
## 
## Number of Fisher Scoring iterations: 5
exp(coef(fit.glm)[1])
## (Intercept) 
##  0.00947787
#confint(fit.glm)
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                           2.5 %     97.5 %
## (Intercept)         0.008143491 0.01101242
## rrh1                0.681558689 0.81398609
## dm1                 1.087462750 1.22331956
## ins_iv1             0.615927037 0.70160385
## ins_sc1             0.644659277 0.73314075
## age                 1.027997904 1.03190001
## gender1             0.828746453 0.92748374
## sofa_1              1.239784611 1.26475456
## elixhauser_hospital 1.082634178 1.09235226
sjp.glm(fit.glm)
## Waiting for profiling to be done...

fit.glm <- glm(mort28 ~ rrh_freq + dm + ins_iv + ins_sc + age + gender + sofa_1 + elixhauser_hospital,
               data=main2, 
               family="binomial")
summary(fit.glm)
## 
## Call:
## glm(formula = mort28 ~ rrh_freq + dm + ins_iv + ins_sc + age + 
##     gender + sofa_1 + elixhauser_hospital, family = "binomial", 
##     data = main2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2450  -0.5555  -0.3882  -0.2445   3.0490  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -4.6635978  0.0769754 -60.586  < 2e-16 ***
## rrh_freq            -0.1557591  0.0291693  -5.340 9.30e-08 ***
## dm1                  0.1428789  0.0300247   4.759 1.95e-06 ***
## ins_iv1             -0.4271564  0.0332210 -12.858  < 2e-16 ***
## ins_sc1             -0.3762760  0.0328080 -11.469  < 2e-16 ***
## age                  0.0295084  0.0009664  30.535  < 2e-16 ***
## gender1             -0.1318217  0.0287112  -4.591 4.41e-06 ***
## sofa_1               0.2236191  0.0050768  44.047  < 2e-16 ***
## elixhauser_hospital  0.0841566  0.0022782  36.940  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 38569  on 48813  degrees of freedom
## Residual deviance: 33312  on 48805  degrees of freedom
##   (9 observations deleted due to missingness)
## AIC: 33330
## 
## Number of Fisher Scoring iterations: 5
exp(coef(fit.glm)[1])
## (Intercept) 
## 0.009432465
#confint(fit.glm)
exp(confint(fit.glm))
## Waiting for profiling to be done...
##                          2.5 %     97.5 %
## (Intercept)         0.00810465 0.01095943
## rrh_freq            0.80763035 0.90548273
## dm1                 1.08759224 1.22344176
## ins_iv1             0.61113702 0.69614164
## ins_sc1             0.64354712 0.73187421
## age                 1.02800568 1.03190763
## gender1             0.82854787 0.92725356
## sofa_1              1.23823245 1.26312255
## elixhauser_hospital 1.08295882 1.09267381
sjp.glm(fit.glm)
## Waiting for profiling to be done...

  1. Stepwise selection
full.model <- glm(mort28 ~ rh + dm + ins + age + gender + sofa_cat + elixhauser_hospital + tpn + vasop + ventilat,
                  data=main2, family="binomial")
stepAIC(full.model, scope=list(lower=~dm)) # AIC
## Start:  AIC=32217.16
## mort28 ~ rh + dm + ins + age + gender + sofa_cat + elixhauser_hospital + 
##     tpn + vasop + ventilat
## 
##                       Df Deviance   AIC
## - tpn                  1    32192 32216
## <none>                      32191 32217
## - gender               1    32210 32234
## - rh                   1    32264 32288
## - ventilat             1    32266 32290
## - ins                  1    32447 32471
## - sofa_cat             3    32677 32697
## - age                  1    33224 33248
## - elixhauser_hospital  1    33433 33457
## - vasop                1    33563 33587
## 
## Step:  AIC=32215.9
## mort28 ~ rh + dm + ins + age + gender + sofa_cat + elixhauser_hospital + 
##     vasop + ventilat
## 
##                       Df Deviance   AIC
## <none>                      32192 32216
## - gender               1    32211 32233
## - ventilat             1    32266 32288
## - rh                   1    32267 32289
## - ins                  1    32448 32470
## - sofa_cat             3    32679 32697
## - age                  1    33233 33255
## - elixhauser_hospital  1    33439 33461
## - vasop                1    33579 33601
## 
## Call:  glm(formula = mort28 ~ rh + dm + ins + age + gender + sofa_cat + 
##     elixhauser_hospital + vasop + ventilat, family = "binomial", 
##     data = main2)
## 
## Coefficients:
##         (Intercept)                  rh1                  dm1  
##            -4.88493             -0.33460              0.10571  
##                ins1                  age              gender1  
##            -0.50553              0.03045             -0.12710  
##      sofa_cat[1, 3)       sofa_cat[3, 6)       sofa_cat[6,15]  
##             0.35352              0.91178              1.07617  
## elixhauser_hospital               vasop1            ventilat1  
##             0.08159              1.89672              0.33804  
## 
## Degrees of Freedom: 48813 Total (i.e. Null);  48802 Residual
##   (9 observations deleted due to missingness)
## Null Deviance:       38570 
## Residual Deviance: 32190     AIC: 32220
stepAIC(full.model, scope=list(lower=~dm), k=log(nrow(main))) # BIC
## Start:  AIC=32331.51
## mort28 ~ rh + dm + ins + age + gender + sofa_cat + elixhauser_hospital + 
##     tpn + vasop + ventilat
## 
##                       Df Deviance   AIC
## - tpn                  1    32192 32321
## <none>                      32191 32332
## - gender               1    32210 32340
## - rh                   1    32264 32394
## - ventilat             1    32266 32395
## - ins                  1    32447 32576
## - sofa_cat             3    32677 32785
## - age                  1    33224 33354
## - elixhauser_hospital  1    33433 33563
## - vasop                1    33563 33693
## 
## Step:  AIC=32321.46
## mort28 ~ rh + dm + ins + age + gender + sofa_cat + elixhauser_hospital + 
##     vasop + ventilat
## 
##                       Df Deviance   AIC
## <none>                      32192 32321
## - gender               1    32211 32329
## - ventilat             1    32266 32385
## - rh                   1    32267 32386
## - ins                  1    32448 32566
## - sofa_cat             3    32679 32776
## - age                  1    33233 33351
## - elixhauser_hospital  1    33439 33558
## - vasop                1    33579 33698
## 
## Call:  glm(formula = mort28 ~ rh + dm + ins + age + gender + sofa_cat + 
##     elixhauser_hospital + vasop + ventilat, family = "binomial", 
##     data = main2)
## 
## Coefficients:
##         (Intercept)                  rh1                  dm1  
##            -4.88493             -0.33460              0.10571  
##                ins1                  age              gender1  
##            -0.50553              0.03045             -0.12710  
##      sofa_cat[1, 3)       sofa_cat[3, 6)       sofa_cat[6,15]  
##             0.35352              0.91178              1.07617  
## elixhauser_hospital               vasop1            ventilat1  
##             0.08159              1.89672              0.33804  
## 
## Degrees of Freedom: 48813 Total (i.e. Null);  48802 Residual
##   (9 observations deleted due to missingness)
## Null Deviance:       38570 
## Residual Deviance: 32190     AIC: 32220

xxx

CreateTableOne(vars=c("ps.groups", "ps.groups2", "age", "gender", "a1c_pre_icu", "sofa_1", "elixhauser_hospital", 
                      "hypoglycemia", "hypoglycemia_freq", "hyperglycemia", "hyperglycemia_freq", 
                      "ins_mean", "tpn", "vasop", "ventilat", "mort28", "mort90", "mort365", 
                      "survival_day", "icu_los", "hosp_los"),
               strata=c("dm"),
               data=main2, test=TRUE) %>% print(
  printToggle      = FALSE,
  showAllLevels    = TRUE,
  cramVars         = "kon"
  ) %>% {data.frame(
  variable_name    = gsub(" ", "&nbsp;", rownames(.), fixed = TRUE), ., 
  row.names        = NULL, 
  check.names      = FALSE, 
  stringsAsFactors = FALSE)} %>% knitr::kable()
## Warning in ModuleReturnVarsExist(vars, data): The data frame does not have:
## ps.groups ps.groups2 Dropped
variable_name level 0 1 p test
n 30755 18068
age (mean (sd)) 62.93 (18.43) 66.09 (14.94) <0.001
gender (%) 0 13408 (43.6) 7957 (44.0) 0.346
1 17347 (56.4) 10111 (56.0)
a1c_pre_icu (mean (sd)) 5.71 (0.41) 7.30 (1.85) <0.001
sofa_1 (mean (sd)) 2.95 (2.83) 2.88 (2.85) 0.014
elixhauser_hospital (mean (sd)) 2.63 (6.27) 2.23 (6.49) <0.001
hypoglycemia (%) 0 26536 (86.3) 12823 (71.0) <0.001
1 4219 (13.7) 5245 (29.0)
hypoglycemia_freq (mean (sd)) 0.28 (1.04) 0.65 (1.57) <0.001
hyperglycemia (%) 0 20125 (65.4) 5035 (27.9) <0.001
1 10630 (34.6) 13033 (72.1)
hyperglycemia_freq (mean (sd)) 0.96 (2.41) 3.77 (5.36) <0.001
ins_mean (mean (sd)) -12.03 (36.22) -5.33 (105.02) <0.001
tpn (%) 0 29168 (94.8) 17225 (95.3) 0.016
1 1587 ( 5.2) 843 ( 4.7)
vasop (%) 0 29397 (95.6) 17109 (94.7) <0.001
1 1358 ( 4.4) 959 ( 5.3)
ventilat (%) 0 15453 (50.2) 9177 (50.8) 0.248
1 15302 (49.8) 8891 (49.2)
mort28 (%) 0 26632 (86.6) 15612 (86.4) 0.568
1 4123 (13.4) 2456 (13.6)
mort90 (%) 0 24984 (81.2) 14407 (79.7) <0.001
1 5771 (18.8) 3661 (20.3)
mort365 (%) 0 22625 (73.6) 12605 (69.8) <0.001
1 8130 (26.4) 5463 (30.2)
survival_day (mean (sd)) 501.58 (751.67) 528.60 (743.80) 0.009
icu_los (mean (sd)) 3.95 (5.83) 4.33 (6.26) <0.001
hosp_los (mean (sd)) 9.78 (10.83) 10.54 (10.83) <0.001